function map = glm_RFX_rMap(hfile, c, r, mapopts)
% GLM::RFX_rMap  - calculate a second-level r contrast map
%
% FORMAT:       map = glm.RFX_rMap(c, r [, mapopts])
%
% Input fields:
%
%       c           NxC contrast vector
%       r           Rx1 or RxN regression data (or VOI object)
%       mapopts     structure with optional fields
%        .const     also create (t-) map of constant term
%        .meanr     boolean flag, remove mean from map (added as cov)
%        .meanrmsk  mask to get mean from (object or XxYxZ logical)
%        .names     1xN cell array with map names
%        .robust    flag, use robust regression in addition to OLS
%        .subsel    subject selection (otherwise all subjects)
%        .thresh    1x2 threshold (lower, upper), as p-values!
%        .voiidx    index into VOI list (only used if r is a VOI object)
%
% Output fields:
%
%       map         MAP/VMP/SMP object with maps

% Version:  v0.7g
% Build:    9092717
% Date:     Sep-27 2009, 5:12 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'glm') || ...
   ~isa(c, 'double') || ...
    isempty(c) || ...
    any(isnan(c(:)) | isinf(c(:))) || ...
   ((~isa(r, 'double') || ...
     isempty(r) || ...
     any(isnan(r(:)) | isinf(r(:))) || ...
    ~any(size(r, 2) == [1, size(c)])) && ...
    (numel(r) ~= 1 || ...
    ~isBVQXfile(r, 'voi')))
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if bc.ProjectTypeRFX ~= 1
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
if ~any(bc.ProjectType == [1, 2])
    error( ...
        'BVQXfile:Unsupported', ...
        'RFX correlation maps of FMRs are not yet supported.' ...
    );
end
numsubs = numel(bc.GLMData.Subject);
if numsubs < 3
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid RFX GLM object.' ...
    );
end
if bc.ProjectType == 1
    msz = size(bc.GLMData.RFXGlobalMap);
else
    msz = numel(bc.GLMData.RFXGlobalMap);
end
numspred = size(bc.GLMData.Subject(1).BetaMaps, ...
    ndims(bc.GLMData.Subject(1).BetaMaps));
if nargin < 4 || ...
   ~isstruct(mapopts) || ...
    numel(mapopts) ~= 1
    mapopts = struct;
end
if ~isfield(mapopts, 'const') || ...
   ~islogical(mapopts.const) || ...
    numel(mapopts.const) ~= 1
    mapopts.const = false;
end
if ~isfield(mapopts, 'meanr') || ...
   ~islogical(mapopts.meanr) || ...
    numel(mapopts.meanr) ~= 1
    mapopts.meanr = false;
end
if isfield(mapopts, 'meanrmsk') && ...
    numel(mapopts.meanrmsk) == 1 && ...
    isBVQXfile(mapopts.meanrmsk, 'msk')
    mbc = bvqxfile_getcont(mapopts.meanrmsk.L);
    if numel(mbc.Mask) == prod(msz)
        mapopts.meanrmsk = lsqueeze(mbc.Mask > 0);
    else
        mapopts.meanrmsk = [];
    end
elseif isfield(mapopts, 'meanrmsk') && ...
    islogical(mapopts.meanrmsk) && ...
    numel(mapopts.meanrmsk) == prod(msz)
    mapopts.meanrmsk = lsqueeze(mapopts.meanrmsk);
else
    mapopts.meanrmsk = [];
end
if isempty(mapopts.meanrmsk) && ...
    mapopts.meanr
    mapopts.meanrmsk = all(bc.GLMData.Subject(1).BetaMaps ~= 0, ...
        ndims(bc.GLMData.Subject(1).BetaMaps));
    for sc = 1:numsubs
        mapopts.meanrmsk = (mapopts.meanrmsk & ...
            all(bc.GLMData.Subject(1).BetaMaps ~= 0, ...
            ndims(bc.GLMData.Subject(1).BetaMaps)));
    end
    mapopts.meanrmsk = lsqueeze(mapopts.meanrmsk);
end
if ~isfield(mapopts, 'names') || ...
   ~iscell(mapopts.names) || ...
    isempty(mapopts.names)
    mapopts.names = {};
end
if ~isfield(mapopts, 'robust') || ...
   ~islogical(mapopts.robust) || ...
    numel(mapopts.robust) ~= 1
    mapopts.robust = false;
end
if ~isfield(mapopts, 'subsel') || ...
   ~isa(mapopts.subsel, 'double') || ...
    isempty(mapopts.subsel) || ...
    any(isinf(mapopts.subsel(:)) | isnan(mapopts.subsel(:))) || ...
    numel(unique(round(mapopts.subsel(:)))) ~= numel(mapopts.subsel) || ...
    any(mapopts.subsel(:) < 1 | mapopts.subsel(:) > numsubs)
    mapopts.subsel = 1:numsubs;
else
    mapopts.subsel = round(mapopts.subsel(:)');
end
if ~isfield(mapopts, 'voiidx') || ...
   ~isa(mapopts.voiidx, 'double') || ...
    numel(mapopts.voiidx) ~= 1 || ...
    isinf(mapopts.voiidx) || ...
    isnan(mapopts.voiidx) || ...
    mapopts.voiidx < 1
    mapopts.voiidx = 1;
else
    mapopts.voiidx = floor(mapopts.voiidx);
end
if numel(r) == 1
    try
        rbc = bvqxfile_getcont(r.L);
        rbc = numel(rbc.VOI);
        if size(c, 2) ~= numspred && ...
            size(c, 2) ~= (numspred - 1)
            ct = c';
        else
            ct = c;
        end
        r = glm_VOIBetas(hfile, r, ...
            struct('c', ct, 'vl', min(mapopts.voiidx, rbc)));
    catch
        rethrow(lasterror);
    end
end
if ~any(size(r) == numel(mapopts.subsel))
    error( ...
        'BVQXfile:BadArgument', ...
        'Correlation regressors must match in size with number of subjects.' ...
    );
end
if ~isfield(mapopts, 'thresh') || ...
   ~isa(mapopts.thresh, 'double') || ...
    numel(mapopts.thresh) ~= 2 || ...
    any(isinf(mapopts.thresh) | isnan(mapopts.thresh) | mapopts.thresh <= 0 | mapopts.thresh >= 0.5)
    mapopts.thresh = [0.005, 0.0001];
else
    mapopts.thresh = -sort(-mapopts.thresh(:)');
end
subsel = mapopts.subsel;
numsubs = numel(subsel);
if mapopts.meanr
    mvm = zeros(numsubs, 1);
end
nval = numsubs - 2;
thresh = mapopts.thresh;
if size(r, 2) == numel(mapopts.subsel) && ...
    size(r, 1) ~= numel(mapopts.subsel)
    r = r';
end
if size(c, 1) == 1 && ...
   (size(c, 2) == (numspred - 1) || ...
    size(c, 2) == numspred)
    c = c';
end
if size(c, 1) == (numspred - 1)
    c(end+1,:) = 0;
end
if size(c, 1) ~= numspred
    error( ...
        'BVQXfile:BadArgument', ...
        'Contrast vector must span all conditions.' ...
    );
end
if size(c, 2) == 1 && ...
    size(r, 2) > 1
    c = c(:, ones(1, size(r, 2)));
elseif size(r, 2) == 1 && ...
    size(c, 2) > 1
    r = r(:, ones(1, size(c, 2)));
end
nummaps = size(c, 2);
nummapst = nummaps * ...
    (1 + mapopts.meanr) * (1 + mapopts.const) * (1 + mapopts.robust);
if bc.ProjectType == 1
    subsa = {':', ':', ':', []};
    subsr = {':', ':', ':', []};
else
    subsa = {':', []};
    subsr = {':', []};
end
rpma = [msz, 1];
if numel(mapopts.names) ~= nummaps
    mapopts.names = cell(1, nummaps);
end
for cc = 1:numel(mapopts.names)
    if ~ischar(mapopts.names{cc})
        mapopts.names{cc} = sprintf('Correlation map %d', cc);
    elseif numel(mapopts.names{cc}) > 96
        mapopts.names{cc} = ...
            [mapopts.names{cc}(1:46) ' ... ' mapopts.names{cc}(end-45:end)];
    else
        mapopts.names{cc} = mapopts.names{cc}(:)';
    end
end

% create map container
switch (bc.ProjectType)
    
    % VTC/VMP
    case {1}
        map = bless(newnatresvmp, 1);
        mapc = bvqxfile_getcont(map.L);
        mapc.XStart = bc.XStart;
        mapc.XEnd = bc.XEnd;
        mapc.YStart = bc.YStart;
        mapc.YEnd = bc.YEnd;
        mapc.ZStart = bc.ZStart;
        mapc.ZEnd = bc.ZEnd;
        mapc.Resolution = bc.Resolution;
        mapf = 'VMPData';
        
    % MTC/SMP
    case {2}
        map = bless(BVQXfile('new:smp'), 1);
        mapc = bvqxfile_getcont(map.L);
        mapc.NrOfVertices = bc.NrOfVertices;
        mapf = 'SMPData';
end

% set some common fields
mapc.Map.NrOfMaps = nummapst;
mapc.Map.Type = 2;
mapc.Map.LowerThreshold = ...
    correlinvtstat(-custom_tinv(thresh(1), nval), numsubs);
mapc.Map.UpperThreshold = ...
    correlinvtstat(-custom_tinv(thresh(2), nval), numsubs);
mapc.Map.DF1 = nval;
mapc.Map.DF2 = 0;
mapc.Map.NrOfFDRThresholds = 0;
mapc.Map.FDRThresholds = zeros(0, 3);
mapc.Map.(mapf) = single(zeros(rpma));

% replicate
mapc.Map = mapc.Map(1, ones(1, nummapst));

% computation
conmaps = zeros([msz, numsubs]);
tmc = 1;
for cc = 1:nummaps
    
    % fill contrast maps
    for pc = 1:numspred
        if c(pc, cc) ~= 0
            subsr{end} = pc;
            subsrs = struct('type', '()', 'subs', {subsr});
            for sc = 1:numsubs
                subsa{end} = sc;
                subsas = struct('type', '()', 'subs', {subsa});
                conmaps = subsasgn(conmaps, subsas, ...
                    subsref(conmaps, subsas) + c(pc, cc) .* ...
                    subsref(bc.GLMData.Subject(subsel(sc)).BetaMaps, subsrs));
            end
        end
    end
    conmaps(isinf(conmaps) | isnan(conmaps)) = 0;
    
    % generate design matrix/ces
    X = [ztrans(r(:, cc)), ones(numsubs, 1)];
    iXX = pinv(X' * X);
    if mapopts.meanr
        for sc = 1:numsubs
            subsa{end} = sc;
            subsas = struct('type', '()', 'subs', {subsa});
            mv = subsref(conmaps, subsas);
            mvm(sc) = mean(mv(mapopts.meanrmsk));
        end
        Xm = [X, ztrans(mvm)];
        iXXm = pinv(Xm' * Xm);
    end
    
    % OLS computations first
    betas = iXX * X' * reshape(conmaps, prod(msz), numsubs)';
    stder = sqrt(sum((conmaps - reshape((X * betas)', [msz, numsubs])) .^ 2, ...
        ndims(conmaps))) .* sqrt(1 / nval);
    
    % first maps
    tmap = betas(1, :)' ./ (stder(:) .* sqrt(iXX(1)));
    tmap(isinf(tmap) | isnan(tmap)) = 0;
    mapc.Map(tmc).(mapf) = reshape(single(correlinvtstat( ...
        tmap, numsubs)), [msz, 1]);
    mapc.Map(tmc).Name = sprintf('%s', mapopts.names{cc});
    tmc = tmc + 1;
    if mapopts.const
        tmap = betas(2, :)' ./ (stder(:) .* sqrt(iXX(2, 2)));
        tmap(isinf(tmap) | isnan(tmap)) = 0;
        mapc.Map(tmc).(mapf) = reshape(single(tmap), [msz, 1]);
        mapc.Map(tmc).Type = 1;
        mapc.Map(tmc).Name = sprintf('%s (intercept-t)', mapopts.names{cc});
        mapc.Map(tmc).LowerThreshold = -custom_tinv(thresh(1), nval);
        mapc.Map(tmc).UpperThreshold = -custom_tinv(thresh(2), nval);
        tmc = tmc + 1;
    end
    
    % with mean removed ?
    if mapopts.meanr
        betas = iXXm * Xm' * reshape(conmaps, prod(msz), numsubs)';
        stder = sqrt(sum((conmaps - reshape((Xm * betas)', [msz, numsubs])) .^ 2, ...
            ndims(conmaps))) .* sqrt(1 / (nval - 1));

        % first maps
        tmap = betas(1, :)' ./ (stder(:) .* sqrt(iXXm(1)));
        tmap(isinf(tmap) | isnan(tmap)) = 0;
        mapc.Map(tmc).(mapf) = reshape(single(correlinvtstat( ...
            tmap, numsubs)), [msz, 1]);
        mapc.Map(tmc).Name = sprintf('%s (mean-rem)', mapopts.names{cc});
        mapc.Map(tmc).DF1 = nval - 1;
        mapc.Map(tmc).LowerThreshold = ...
            correlinvtstat(-custom_tinv(thresh(1), nval - 1), nval + 1);
        mapc.Map(tmc).UpperThreshold = ...
            correlinvtstat(-custom_tinv(thresh(2), nval - 1), nval + 1);
        tmc = tmc + 1;
        if mapopts.const
            tmap = betas(2, :)' ./ (stder(:) .* sqrt(iXXm(2, 2)));
            tmap(isinf(tmap) | isnan(tmap)) = 0;
            mapc.Map(tmc).(mapf) = reshape(single(tmap), [msz, 1]);
            mapc.Map(tmc).Type = 1;
            mapc.Map(tmc).DF1 = nval - 1;
            mapc.Map(tmc).Name = sprintf('%s (intercept-t, mean-rem)', mapopts.names{cc});
            mapc.Map(tmc).LowerThreshold = -custom_tinv(thresh(1), nval - 1);
            mapc.Map(tmc).UpperThreshold = -custom_tinv(thresh(2), nval - 1);
            tmc = tmc + 1;
        end
    end
    
    % compute robust stats
    if mapopts.robust
        
        % perform fit
        [b, w] = fitrobustbisquare_img(X, conmaps);
        rt = robustt(X, conmaps, b, w, [1, 0]);
        rt(isinf(rt) | isnan(rt)) = 0;
        corm = correlinvtstat(rt, numsubs);
        corm(isinf(corm) | isnan(corm)) = 0;
        mapc.Map(tmc).(mapf) = single(corm);
        mapc.Map(tmc).Name = sprintf('%s (robust)', mapopts.names{cc});
        tmc = tmc + 1;
        if mapopts.const
            rt = robustt(X, conmaps, b, w, [0, 1]);
            rt(isinf(rt) | isnan(rt)) = 0;
            mapc.Map(tmc).Type = 1;
            mapc.Map(tmc).(mapf) = single(rt);
            mapc.Map(tmc).Name = sprintf('%s (robust, intercept-t)', mapopts.names{cc});
            mapc.Map(tmc).LowerThreshold = -custom_tinv(thresh(1), nval);
            mapc.Map(tmc).UpperThreshold = -custom_tinv(thresh(2), nval);
            tmc = tmc + 1;
        end
        if mapopts.meanr
            [b, w] = fitrobustbisquare_img(Xm, conmaps);
            rt = robustt(Xm, conmaps, b, w, [1, 0, 0]);
            rt(isinf(rt) | isnan(rt)) = 0;
            corm = correlinvtstat(rt, numsubs);
            corm(isinf(corm) | isnan(corm)) = 0;
            mapc.Map(tmc).(mapf) = single(corm);
            mapc.Map(tmc).Name = sprintf('%s (robust, mean-rem)', mapopts.names{cc});
            mapc.Map(tmc).DF1 = nval - 1;
            mapc.Map(tmc).LowerThreshold = ...
                correlinvtstat(-custom_tinv(thresh(1), nval - 1), nval + 1);
            mapc.Map(tmc).UpperThreshold = ...
                correlinvtstat(-custom_tinv(thresh(2), nval - 1), nval + 1);
            tmc = tmc + 1;
            if mapopts.const
                rt = robustt(Xm, conmaps, b, w, [0, 1, 0]);
                rt(isinf(rt) | isnan(rt)) = 0;
                mapc.Map(tmc).Type = 1;
                mapc.Map(tmc).(mapf) = single(rt);
                mapc.Map(tmc).Name = sprintf('%s (robust, intercept-t, mean-rem)', mapopts.names{cc});
                mapc.Map(tmc).DF1 = nval - 1;
                mapc.Map(tmc).LowerThreshold = ...
                    -custom_tinv(thresh(1), nval - 1);
                mapc.Map(tmc).UpperThreshold = ...
                    -custom_tinv(thresh(2), nval - 1);
                tmc = tmc + 1;
            end
        end
    end
    
    % zero out conmaps if necessary
    if cc < nummaps
        conmaps(:) = 0;
    end
end

% put back
bvqxfile_setcont(map.L, mapc);
